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ABSTRACT 


This research was conducted to enhance understanding of the use of 
high explosives to simulate the effects of a nuclear underwater explosion. A 
review of the known characteristics of the nuclear, spherical conventional, 
and tapered conventional underwater pressure-time histories illustrates the 
selection of the tapered charge to simulate the underwater nuclear explosion. 
Three areas of study were then pursued. The first compared the structural 
response resulting from attack by conventional and nuclear type pressure 
profiles, verifying the need to match duration as well as peak pressure when 
simulating the underwater nuclear explosion. The second employed finite 
element analysis to study the three dimensional shock generated by a 
tapered charge. Third, a computer program was written to couple an 
optimizer with an existing tapered charge pressure-profile generating code to 


improve the tapered charge design process. 
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I. INTRODUCTION 


This thesis addresses efforts to improve simulation, using conventional 
high explosives and scale models, of the underwater shock environment and 
structural effects resulting from the underwater nuclear explosion. 
Simulation using scale models and small conventional charges provides 
valuable information without the requirement for nuclear testing, with 
minimal environmental impact, and at low cost. Better understanding of the 
physics involved impacts ship and weapons designs. 

In the current atmosphere of reduced military spending in order to reap 
the benefits of the "Peace Dividend" resulting from the end of the cold war, 
the threat to U. S. Navy ships and submarines from underwater nuclear 
explosion would appear to be greatly reduced. However, two factors ensure 
the continued existence of the threat from underwater nuclear explosion: 

1. The presence of the former Soviet Union's vast arsenal of nuclear 
weapons combines with economic instability to increase the likelihood of 


more nations gaining access to the material necessary to construct 
nuclear weapons. 


2. The ceaseless march of technology worldwide dictates future 
growth in the number of nations attaining the particular technology 
necessary to build and detonate nuclear weapons. 
The growing community of nuclear weapons capable nations may not posses 
the same restraint from the use of nuclear weapons displayed since 1945. 
Given the existence of a threat from underwater nuclear explosion, the 
significance of this threat can be determined only if the effects are well 


understood. This same understanding is essential to incorporation of shock 


hardening in ship and submarine designs to improve survivability. 


With the overall objective of improving ship and submarine 
survivability through better understanding of the phenomena associated 
with underwater explosions, the Naval Postgraduate School conducts 
ongoing research into underwater explosions and effects. This thesis is the 
result of a part of that continuing research, the first at the school related 
specifically to the nuclear underwater explosion. 

The known characteristics of the nuclear underwater explosion, 
together with a discussion of modeling techniques is found in Chapter II. 

Chapter III presents a comparison of the structural response of a simple 
cylinder subjected to side-on attack by conventional spherical charge and 
nuclear type pressure profiles. The doubly asymptotic boundary assumption 
combined with an explicit finite element method was used to perform the 
analysis. Results verify the need to match peak pressure and duration of the 
nuclear pressure-time history when designing test charges to simulate the 
underwater nuclear explosion. 

The tapered charge, due to inherent long pressure duration, commonly 
generates the simulated nuclear pressure field used in model testing. 
Chapter IV explores the three dimensional aspects of the shock front 
generated by a conventional underwater tapered charge explosion. The 
results of analysis performed using an explicit finite element method are 
presented with the intent of supplementing existing knowledge of the 
tapered charge pressure profile which to date consists mostly of on-axis data. 
The results show the evolving shape of the shock front at early times and 
provide information on the relationship between the peak pressures 


measured on and off the charge axis. 


As outlined in Chapter V, a computer code was wnitten to optimize the 
design process for a conventional tapered high explosive charge. Starting 
with a desired pressure-time history and initial estimates for charge 
geometry and standoff distance, this program utilizes public domain 
optimization software to return improved design values. Although tested 
with a subroutine based upon an existing routine for calculating the 
pressure-time history of a tapered charge, this program may be coupled with 
other existing or future routines for calculating the tapered charge 
pressure-time history. 

Conclusions and recommendations for further research 1n this area may 


be found in Chapter VI. 


Il. CHARACTERISTICS AND MODELING OF THE 
UNDERWATER NUCLEAR EXPLOSION 


Before attempting to simulate the underwater nuclear explosion using 
conventional high explosive charges for scale model testing, a degree of 
familiarity with the shock generated by nuclear and spherical high explosive 
charges 1s warranted. This chapter, therefore, outlines the use of empirical 
relationships to determine the pressure-time histories generated by the two 
types of charges. The tremendous weight and standoff distance required for 
a conventional spherical charge to create a pressure profile similar to that of 
a nuclear charge points directly to the need for scaling. 

Following a brief description of the principles used to construct scale 
models for underwater shock testing, these principles are applied to nuclear 
and conventional spherical charges. The resultant large size and standoff, 
even after scaling, of the conventional spherical eres leads to the selection 
of the tapered charge, made of conventional high explosives, to simulate the 


underwater shock generated by the underwater nuclear explosion. 


A. THE UNDERWATER NUCLEAR SHOCK PROFILE 

The energy content, or "yield", of a nuclear explosion is commonly 
measured in tons of TNT equivalent, the amount of explosive energy 
contained in 2,000 pounds of the conventional high explosive TNT. A one 
kiloton (kT) device contains the energy equivalent of 1,000 tons of TNT, one 
megaton (MT) that of 1,000,000 tons. Although the majority of the energy 


released in an underwater nuclear explosion contributes to the generation of 


the underwater shock wave, the extremely high temperatures (tens of 
millions of degrees) reached in a nuclear explosion contribute to a significant 
amount of energy release in the form of thermal radiation. Chemical 
explosions, by contrast, occur at much lower temperatures (thousands of 
degrees), resulting in a higher percentage of the total energy released as 
kinetic energy to generate the underwater shock. (Glasstone and Dolan, 
1977, pp. 1-3, 6, 11) 

To date, the United States has conducted five announced underwater 
nuclear explosions, from 1946 to 1962 (Bolt, 1976, pp. 251-274). Glasstone 
and Dolan (1977, pp. 268-272) provide three empirical charts to calculate the 
pressure-time history of an underwater nuclear explosion given the yield of 
the device and the standoff distance R from the explosion. From the curves 


of the first two charts, the maximum pressure P | (psi) and the time 


constant 6 (ms) are determined. The time constant equals, as in exponential 
decay, the time between the arrival of the shock when P = P_.. and the time 
at which P = P_ /e = 0.37P_... The pressure actually decays at a somewhat 
slower than exponential rate after one time constant, necessitating the third 
chart which plots the non-dimensional values P(t)/P_.. vs t/9 to provide an 
idealized pressure-time history for an underwater nuclear explosion with no 
bottom or surface reflection effects. Based upon these curves, Figure 1 shows 
the pressure time history of a 40 kT nuclear explosion at a standoff of 1,000 
yds. Figure 2 illustrates standoff. Two prominent features of the 


underwater nuclear explosion stand out in Figure 1: the high pressure at a 


significant distance from the explosion, and the long decay time. P... 


increases with yield, decreases with standoff: the duration of the shock wave 


increases with yield and standoff (Glasstone and Dolan, 1977, p. 269). 




















Figure 1. Pressure profile of a 40 RT nuclear charge, R — 1000 yds. 
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Figure 2. Measurement of standoff distance R. 


In water of sufficient depth where bottom reflections are negligible or 
occur at a much later time, the steady pressure decay ends abruptly by the 
phenomenon of surface cutoff (Shin and Geers, 1991, $3.3; or Glasstone and 
Dolan, 1977, pp. 244-246). Figure 3 shows two paths followed by a shock 
wave emanating from an underwater explosion. The direct, compression 
wave strikes the target first with a sudden rise in pressure followed by the 


steady pressure decay described above. The other path shows a compression 


wave striking the air-water interface and being reflected as a tension, or 


rarefaction wave. 


Air 
Water 





Compressive Shock Wave > 


Charge Target 
Figure 3. Paths followed by the direct and rarefaction waves. 


Approximating shock speed by the speed of sound in water C, gives the 
cutoff time t. as the difference in distance traveled by the direct and 


rarefaction waves divided by C: 


2 |(Е) «D? -R 


tco — 
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The arrival of the rarefaction wave at time t,, after the arrival of the 
direct shock wave causes a sudden pressure drop as the remaining pressure 


from the compression wave is essentially canceled. Hence, if the 40 kT 


~] 


nuclear charge discussed earlier is detonated at a 285 ft depth in deep water, 
a target at 1,000 vards and the same depth will experience t,, = 10.7 ms using 
the expression on the previous page with C, — 5 ft/ms. Тһе resultant 


pressure profile with pressure cutoff is shown in Figure 4. 

















Figure 4. Pressure profile of a 40 RT nuclear charge with surface cutoff. 
R=1,000 yds, D=285 ft. 


For an actual underwater nuclear explosion, the smooth curve of Figure 
4 would be further modified bv refraction of the shock wave due to salinity, 
currents, and temperature variation in the water media. Whereas Figure 4 
depicts an instantaneous pressure rise upon arrival of the compression wave 
at time 0, a finite rise time would be expected. Additionally, the abrupt 
surface cutoff shown assumes the rarefaction wave and the compression 
wave travel at the same speed of sound in water. In actuality, the 
rarefaction wave travelling in shocked media partially overtakes the 
compression wave rendering a less steep pressure drop at t. For explosions 
occurring in shallow water, bottom reflections and retransmissions further 


alter the shock profile. Additional shocks may occur at later times due to the 


bubble pulse resulting from explosions at depths such that the gases of the 
explosion expand and collapse before venting at the surface. More than three 
bubble pulses are unlikely due to steam condensation. (Glasstone and Dolan, 
1977, pp. 56, 245, 246, 269) 

Having established the general characteristics of the pressure-time 
history of the underwater nuclear explosion, the next section uses empirical 
formulas to determine the feasibility of emulating this profile using 


conventional high explosives of spherical shape. 


B. SPHERICAL CONVENTIONAL CHARGE PROFILE 

As in the case of the nuclear charge, the conventional charge of 
spherical or near spherical shape gives rise to a sudden pressure increase 
followed by exponential pressure decay for one time constant 0 and somewhat 
slower decay thereafter. Using exponential pressure decay for an 
approximation, empirical studies provide the following useful relationships 
for determining the pressure-time history developed by a conventional high 
explosive of spherical or near spherical shape when detonated underwater 


(Shin and Geers, 1991, §3.2): 





P(t) = | 
_ (па | | 
а= к( R ) = maximum pressure (psi) 





1/3 Аз | De. 
Сее e = time (ms) for P to decay to —— 


where K;, Ko, Aj, A» - empirical constants for a given explosive 
t = time (ms) 
P = pressure (psi) 
R = standoff (ft) 
W = charge weight (lb). 


Using the empirical relationships above, a charge of 27.5 thousand tons 
of TNT would be required to emulate the pressure- time history generated by 


a 40 kT underwater nuclear explosion as illustrated in Figure 5. 


























Figure 5. Comparison of 40 kT nuclear (solid line) and 27,500 ton TNT 
(dashed line) explosions. R=3,000 ft, D=285 ft for both charges. 

The fact that 27.5 kT of TNT generates a similar pressure profile as a 
nuclear device of 40 kT yield is attributed to the thermal energy released by 
the nuclear explosion as discussed in Section II.A. Obviously, using 27.5 
thousand tons of TNT to study the effects of a nuclear explosion is not 
practical, nor feasible. This amount of high explosive represents a volume 
greater than that of the Washington Monument (Lapp, 1980). Scaling laws, 
discussed in the next section, enable the use of scale models and charges of a 
more reasonable size to simulate the underwater nuclear explosion and its 


effects. 
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C. HOPKINSON SCALING 

Dimensional analysis yields scaling principles which enabie the use of 
scale models to replicate the behavior of full size objects, or prototypes. Of 
particular utility in the analysis of underwater explosions and effects is 
Hopkinson scaling. Through the use a scale factor à, the quantities length, 


mass, and time are scaled as follows: 


Model Prototype 
AL Length L 
At Time t 
AM Mass M 


giving the invariant quantities: 


p Density p 
с, Р Stress, Pressure с, Р 
Є Strain e 


Some quantities, chiefly the hydrostatic loading due to gravity, are not 
adaptable to scaling and require additional consideration in modeling. (Shin 
and Geers, 1991, §4.1) 

The benefits of scaling are many. By reducing standoff and charge size, 
tests can be conducted in small manmade ponds with little or no 
environmental impact. Geometrically similar models can be constructed of 
the same materials as the prototype to study structural response to 
underwater shock. Model stress and strain levels will match those of the 
prototype. Since the material required to build model charges and structures 
is equal to à” that of the prototype, model testing delivers obvious cost 


benefits. 


ІШ 


Using a scale factor of A = 1/30 to simulate the 40 kT nuclear 
pressure-time history discussed previously would require a spherical TNT 


charge of size 


Ibm 


ton 


9 
Wry =A°Wpe= E) x 28.5 x 10° x tons x2,000 


The standoff required is 


Кый = T x 3,000 ft — 100 ft. 


Figure 6 shows the 40 kT nuclear pressure-time history scaled using 4 = 1/30 


= 2,040 lbm. 


together with the pressure profile of 2,040 pounds of TNT. The only 


difference between this figure and Figure 5 is the time scale. 











Figure 6. 2,040 lbm TNT charge (dashed line), R=100 ft, to simulate a 
40 kT nuclear (solid line) explosion using a scale factor of A = 1/30. 


Although the spherical conventional charge size is now feasible, the 
weight and standoff required for the simulation are too great for small pond 
testing limited to nominal charge weights and standoffs in the tens of pounds 


and feet respectively. One might be tempted in scale model testing to use a 
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smaller charge and shorter standoff to match only the peak pressure of the 
nuclear profile when studying shock effects on a scale model. As shown in 
Figure 7, a 14 pound spherical TNT charge at a standoff of 19 ft gives the 
same maximum pressure as a 2040 pound charge at a standoff of 100 ft. The 
pressure of the lighter charge decays much more rapidly due to its smaller 


time constant. 





Figure 7. Pressure-time histories of a 2,040 lbm (solid line) TNT charge 
at R=100 ft with a 14 lbm (dashed line) TNT charge at R=19 ft. 

It will be shown in Chapter III that the structural response from two 
pressure profiles having the same peak pressure but different duration of 
high pressure give rise to dramatically different structural responses. 
Therefore, in order to simulate the slow pressure decay of the underwater 
nuclear explosion using conventional charges of modest size, the shape of the 
conventional charge must be modified. The resulting shape is that of the 


tapered charge discussed in the next section. 
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D. TAPED CHARGE PRESSURE PROFILE 

Designed to generate long duration shock waves to simulate nuclear 
underwater shock loading on scale models, tapered charges consist of a series 
of truncated cones on a common axis fitted with a detonator on the small or 
nose end as shown in Figure 8. Constructed in sizes ranging from a few 
ounces to over 15,000 pounds, the tapered charge generates a directional 
pressure field with maximum duration along the nose side on the charge 


axis. (Gordon and Davidson, 1983) 


tail detonator 








Figure 8. Geometry of the tapered charge. 


For tapered and spherical charges of the same weight and at the same 
standoff, Figure 9 illustrates the trade-off between peak pressure and 
duration distinguishing the two designs. The tapered charge creates a peak 
pressure at a lower value than the spherical charge, followed by a more 


gradual pressure decline over a region called the pressure plateau. At the 
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end of the pressure plateau, the pressure-time history of the tapered charge 


gives way to exponential decay then surface cutoff. 


3500 
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Figure 9. Pressure profiles of tapered (solid line) and spherical (dashed 
line) charges with the same weight and standoff. 

The pressure plateau generated by the tapered charge gives an obvious 
advantage over the spherical charge in simulating the nuclear shock profile. 
The duration of the pressure plateau t, can be determined to a first 
approximation by estimating the time difference between (1) travel from nose 
to tail along the charge axis a distance L at the detonation speed C, then a 
distance L + R to the target at the speed of sound in water C, , and (2) travel 
from the nose to the target a distance R at speed C, . The resulting 


expression for pressure plateau duration is 


_{L ,L+R}| R_,;f1,1) 
"P eR И | 





(РЕН С. 
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The above formula predicts a plateau duration independent of standoff. As 
Figure 10 shows, however, the duration of the pressure plateau actually 


decreases somewhat with increasing standoff. 
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Figure 10. Tapered charge pressure profiles at varying standoff. 


Figure 10 also shows the approximate 1/R relationship between peak 
pressure and standoff. In addition to varying the standoff to match peak 
pressures and the overall length to match pressure plateau duration, the 
segment lengths and diameters of the tapered charge can be varied to taylor 
the shape of the tapered charge pressure-time history to model a particular 
nuclear pressure profile. Design of a tapered charge involves first examining 
existing data for a charge design which produces a pressure-time history 
most nearly matching that desired. The charge design is then adjusted using 


computer calculations tempered with the experience of the designer. Small 
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scale experiments may be used to verify the design before production of a 
large tapered charge. (Costanzo, 1991) 

Additional knowledge of the tapered explosion as well as any 
streamlining of the design process should improve the design of tapered 
charges used to simulate the underwater nuclear explosion. An application 
of the finite element method (FEM) to study the early time propagation of the 
shock generated by a tapered charge explosion is found in Chapter IV. 
Chapter V outlines the application of computer design optimization 


techniques to enhance the tapered charge design process. 
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III. EFFECT OF PRESSURE DURATION ON 
STRUCTURAL RESPONSE 


As mentioned in Section II.C, if one could match only the peak pressure 
when conducting a simulation of the underwater nuclear attack, a spherical 
charge of modest size would suffice. As illustrated in this chapter, however, 
matching of peak pressures alone is not adequate. Comparison of the 
response of a simple cylindrical shell to side-on attack by a long duration, 
nuclear type, and a short duration, conventional type pressure profile yielded 


substantially different results using numerical techniques. 


A. ATTACK CURVES 

In order for a realistic comparison in the model testing environment, the 
pressure profiles used in this study were generated from two 56 lbm HBX-1 
charges. The short duration pressure profile was derived from the empirical 
relationships of Section II.B for a spherical charge at a standoff of 20 ft. The 
long duration pressure profile used in this study to simulate the nuclear 
profile was derived from scaled tapered charge data. Figure 11 shows the 
two attack curves with the same peak pressure of approximately 3,400 psi 
but very different high pressure duration times. Each pressure-time history 
with corresponding standoff was entered into an existing computer code to 
provide underwater shock loading of a simple cylindrical aluminum shell 


model as described in the next section. 
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Figure 11. Two attack curves used for comparison of structural 
response. 


B. NUMERICAL MODEL 

The attack geometry for this study 1s shown 1n Figure 12. The figure 
illustrates the standoff for a spherical charge. Tapered charge standoff 1s, by 
convention, measured from the nose, Figure 8, to the target. The pressure 
profiles and standoffs from the previous section were used with a side-on 
attack cylindrica] shell model developed by Fox (1992. pp. 60.61). The 
analysis was conducted using the public domain finite element method 
(FEM) code VEC/DYNA3D (Hallquist and Stillman, June, 1990) coupled with 
the boundary element method code USA (Deruntz, 1989). The USA 
(Underwater Shock Analyzer) code reduces the media surrounding the 
cylinder and the associated forces to discrete forces and masses to provide 
loading to the cylinder. The FEM code VEC/DYNASD utilizes explicit time 


integration to provide the response of the cylinder to the applied shock 
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loading. The coupling of the two codes was initiated by the Naval 
Postgraduate School (Fox, Kwon, and Shin, 1992). 


42" 


Target Charge 


– 1» |- 


Figure 12. Attack geometry for cylindrical shell FEM experiment. 


The target consisted of a 1/4" thick 6061-T6 aluminum cylinder with 1" 
machined and welded endplates. Material properties for the cylinder were: 


o, = 40,000 psi 
E = 10,800 kpsi 
у = 0.33 

р = 174 lbm/ft’. 

Figure 13 depicts the 550 shell elements comprising the one-quarter 
symmetry discretization of the cylindrical shell model used in this study. A 
description of the theory behind the shell elements used can be found in the 
article by Belytschko, Lin, and Tsay (1984). Symmetric boundary conditions 
were applied on the yz and zx planes. The boundary element loading 
described previously provided the boundary conditions for the outer surface 


of the wet elements of the aluminum cylinder. 
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element 211 


element 214 


Figure 13. Quarter symmetry FEM model of simple cylindrical shell. 


The aluminum was modeled as a kinematic/isotropic/elastic/plastic. In 
this 1dealization. elastic deformation occurs at stress levels lower than the 
yield stress with plastic deformation occurring at the yield stress. This 
provides a reasonable approximation of the uniaxial stress-strain curve for 
6061-T6 aluminum based upon pull-test data as shown in Meyers and Murr 
(1981, p. 40), with the exception of no provision for failure once a maximum 
engineering strain of 7 to 9 per cent has been sustained. Taking 8 per cent 
plastic strain as a failure criterion and applying a factor of safety of 2, failure 
of the aluminum shell model was predicted for effective plastic strain in 


excess of 4 per cent Effective plastic strain £, is defined as (Ugural and 
Fenster, 1987): 


Ae e ea] 
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where ¢,, €., and е, аге the true plastic strain components. Since the 
maximum strain-rate calculated in this study was approximately 100 
in/in-sec, well below the 2000 in/in-sec required for appreciable strain-rate 
effects, (Meyers and Murr. 1981, p. 50), strain-rate hardening was not 


included in the model. 


С. RESULTS 

Figures 14 and 15 show the damage sustained by the cylindrical target 
as a result of attack by the nuclear and conventional pressure profiles with 
the same peak pressures. The quarter-symmetry model has been reflected 
about the yz and zx planes using the post-processor TAURUS (Hallquist, 
1990) to provide visualization of the full model in each case. Displacements 
for the conventional attack of Figure 14 were scaled by a factor of 10 to better 
display the damage pattern. No scaling was done in Figure 15 where the 


damage was much more extensive due to the nuclear type attack. 
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Figure 14. Effective plastic strain from conventional attack (displace- 
ments scaled by a factor of 10). 
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Figure 15. Effective plastic strain from nuclear attack. (No scaling of 
displacements) 


Figure 16 illustrates the time history of effective plastic strain for the 
elements sustaining maximum damage from each attack. The maximum 
damage of 1.67 per cent effective plastic strain experienced by the 
conventional attack, although significant, did not exceed the failure criteria 
selected in the previous section. The effective plastic strain from the nuclear 
type attack, exceeding 4 per cent after less than 0.6 ms, indicates the 
prediction of catastrophic failure of the cylinder resulting from this attack 
based upon the numerical analysis. 

Not only is the magnitude of effective plastic strain much greater for the 
nuclear type attack, the mode of damage is quite different. In both cases. 
maximum damage occurred in elements located approximately 3.2 inches 
from the endplate. In the conventional case, the maximum damage occurred 
in element 211, on the yz plane. Maximum damage in the nuclear case, 


however, occurred at element 214, 30 degrees off the yz plane. 
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however, occurred at element 214, 30 degrees off the yz plane. The locations 


of the two critical elements are shown on Figures 13 through 15. 
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Figure 16. Plot of effective plastic strain for the elements sustaining 
maximum damage. 


Based upon the analysis of this chapter, both the peak pressure and 
high pressure duration must be generated by the бб она! high explosive 
used to simulate the underwater nuclear explosion. Due to the ability of the 
tapered charge to accomplish this task at a lower charge weight than the 
spherical charge, the studies of Chapters IV and V were conducted to 


enhance understanding and improve the tapered charge design process. 
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IV. THE TAPERED CHARGE SHOCK FRONT 


Due to its ability to develop a long duration pressure profile from a 
modest charge weight, the tapered charge commonly provides the shock 
loading to simulate the underwater nuclear attack as discussed previously. 
In order to better understand the shock developed by tapered charges, finite 
element analysis (FEA) was conducted to study the shock developed by a 
tapered charge detonated underwater. The FEA provided information on the 
directional nature of the pressure-time history generated by the tapered 
charge. This information, specific for the charge geometry and type of 
explosive, can be used to determine the accuracy required to position the 


charge and target used in model studies of the underwater nuclear explosion. 


A. FEA MODEL 

Figure 17 shows the 46" x 46" x 152" quarter-symmetry FEA model used 
for this study. The mesh consists of 120 HBX-1 charge elements and 40.692 
water elements. Appendix A provides a listing of the input file for the 
INGRID (Stillman and Hallquist, 1991) mesh generating program used. The 


dimensions of the tapered charge, corresponding to those of Figure 8, were: 


ИІ 00335 $6 


es = 4.330 it 
L = 5.000 ft 
di=1125in 
do = 27625 in 
da — 4.125 in 
G@f= 9010 1. 
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Using a cast specific gravity for HBX-1 of 1.71 (Dobratz, 1981. p. 19.53) gives 
the charge weight of 60.1 lbm. Referring to Figure 17, symmetric boundary 
conditions were applied on the yz and zx planes. Non-reflective boundary 


conditions were applied to the remaining four planes. 
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Figure 17. FEA model used to examine the three dimensional aspects of 
the tapered charge. 

The charge and water model were analyzed using FEA program 
VEC/DYNA3D (Hallquist and Stillman. 1990). This is the stand-alone 
version of the coupled program used in Chapter III. The explosive was 
modeled as a high explosive burn material with the Jones-Wilks-Lee (JWL) 


equation of state, the water as a null material using the Gruneisen equation 


of state. 
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For the charge, the high explosive burn material model used by 
VEC/DYNAS3D requires entry of the detonation velocity D, the Chapman- 
Jouget pressure P..,, and the density p (Hallquist and Stillman, 1990, p. 40). 
The values О = 0.731 cm/s, P., = 0.2204 Mbar, and p = 1.712 gm/cm’ for 
HBX-1 were taken from Dobratz (1981, p. 19.53). The JWL equation of state 
was used to describe the pressure-volume-energy behavior of the detonation 


products (Dobratz, 1981, p. 8.21): 


В, Е ey B -Rov | OF 
ОЕ А(1 R, mos +Ва RE - V 
Where A, B, and C = linear coefficients in Mbar 


Rı. R2, ando = nonlinear coefficients 
volume of detonation products 
А = = ит ле ис е т ыз лыла ж s 
| volume of undetonated high explosive 
P = pressure in Mbar | 
Маг х ст” 


E = det onation energy per unit volume in oma 


The parameters E, A, B, R.. R., and o listed above, empirically derived from 
cylinder-test data, are required entries for use of the JWL equation of state 
with VEC/DYNA3D (Hallquist and Stillman, 1990, p. 89). For this study 
these values were taken from cylinder-test data for H-6, an explosive of 
similar composition to HBX-1, due to non-availability of HBX-1 data. The 
values used were (Dobratz, 1981, p. 8.22): 


3 
= og parcem” 
ст: 


A= 7.5807 Mbar 
B = 0.08513 Mbar 


ЕК, =4.9 
КӨШІ 
wo = 0.20. 


For the water, the VEC/DYNA3D null material model (Hallquist and 
Stillman, 1990, p. 41) requires entry of the density and an optional pressure 
cutoff. The value p - 1.000 gm/cm? was used for density, and a pressure 
cutoff value of 6.89 x 10? Mbar (0.1 psi) were used since water is unable to 
sustain tension. The Gruneisen equation of state with cubic shock velocity 
-particle velocity (u,-u,) defines pressure p in Kbar for compressed materials 
as (Hallquist and Stillman, 1990, p. 91): 

оС? 1 + (1 - a TES su | 
p=: ~ аи (yo + au) E 
u2 TES 
1- (S; - Jusis 
p 


where и soe 1 


. . gm 
= density in 
p nsity E 





. • 1 
Oo = standard density in n 
cm: 


and the required parameters are 


С = е intercept of the u.-u, curve in jus) 
S1, S», Sa- coefficients of the slope of the u.-u, curve 
уо = the Gruneisen gamma 
a = first order volume correction to yo. 
3 
E = internal energy per unit volume in Mar cem. 


The parameters C, S,, S,, and S, for the u,-u, curve were obtained from 
Steinberg (1987): 
Uy uy? 
us =C+S uy disini + ва (12) Up 


- 0.148 4 2.56u, - 1.986, ^u, + 0.2268( <= ) “Up 
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Appendix B outlines the determination of y, = 0.4934 and a = 1.3937, based 
upon the seventh order polynomial approximation for the Gruneisen gamma 


as a function of specific volume developed by Gurtman, Kirsch, and Hastings 


(1971). 


B. RESULTS 

The large number of elements led to a computationally intensive finite 
element analysis of the tapered charge underwater explosion. Running on a 
UNIX engineering workstation required approximately four days to perform 
the computations through 2 us after detonation. Figure 18 shows 
representative time histories of two elements located at a standoff of 35" from 
the charge nose. One element is on the charge axis off the nose, the other 90 
degrees off this axis. Due to the large variation in magnitude of pressures at 
the two locations. as will be shown in Section IV.B.2. the pressures have been 
normalized to better compare the general shapes of the curves. 

The element located on the charge axis maintained a greater portion of 
its maximum peak pressure over a longer period of time than did the element 
located 90 degrees off the axis. This observation further supports the use as 
well as the orientation of the tapered charge for simulation of the underwater 
nuclear explosion. 

Both of the curves of Figure 18 rose less rapidly than expected and 
displayed oscillatory behavior. Gordon and Davidson (1983) experienced 
similar results when analyzing a tapered pentolite charge using 
finite-difference techniques in two dimensions. They attributed the long rise 
time to three possible causes: the artificial viscosity coefficient built into 


their model, the mesh size. and the fact that their model was approximated 
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by a pointed nose. Although the model used in this study did not use a 
pointed nose approximation, the other two possible causes apply to this 
study. Fox (1992. pp. 9.10) cited mesh reflection effects as a contributing 
factor to instabilities when using the FEA code of this study. 








| On Charge Axis | 
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US | | 





























Figure 18. Normalized Pressure-time History of two elements located at 
R= 45 inches: on the charge axis and 90 degrees off the charge axis. 


Although desirable, a finer mesh was beyond the capabilities of the 
machine used in this analysis. For this study, the default artificial viscosity 
coefficient and time integration steps (Hallquist and Stillman, 1990, p. 25) 
were used. Future sensitivity studies may identify values for these 
parameters. The mesh reflection effect appears to be a strong factor in the 
observed pressure oscillations. The key contributor to this problem, uneven 
size of adjacent elements, was aggravated in this study by the necessity to 


conform the shape of eight-node elements to the curvatures of the tapered 
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charge. The end result was fine charge elements near the nose next to 
relatively coarse water elements and a reversal of this effect at the tail of the 
charge. 

Oscillations in pressure and the limited time of the pressure histories 
generated precluded determination of the pressure plateau duration. The 
remainder of this section concentrates on the shock wave as it emanates from 
the charge at early times and the directional nature of the peak pressures in 
the media surrounding the tapered charge. 

1. Early Time Shape of the Shock Front 

Figures 19 through 25 show a view perpendicular to the yz-plane 
of the water elements. The model has been reflected about the zx-plane, with 
charge elements removed. Figures 19 applies to time pror to detonation. 


Shading indicates the elements used for the plots of Section IV.B.2. 





Figure 19. Water elements prior to detonation. 
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Figures 20 through 25 comprise a set of pressure contour plots at 
early times through 0.58 ms. There are five contour lines in each plot. The 
contours range from 1,000 psi for the sparsest dotted line to 5,000 psi for the 
solid line. This relatively low pressure range at such close proximity to the 


charge was selected to clearly define the location of the evolving shock. 






EEEE I BET 
Рг aeaaeae ЕЕ оиро 
ш и АДАД ДЕДА ДАДЕ ЕАУ АЕ РРА РАДАР ааа и ои“ 
ДАДА АДАД ЕРЕ ааг а аашаа ааа || [|] ее ш ааа НЕ 
ШӘЕРЕГІГІІШГІІТЕРІТІРГРІЕТГТІЕРІГІТІТІГІГІТІРІГІРІІГІЕРГРІТТІГІШРІРІІТТІГІТГІІЕГРІ раан ие 
SS Se SES уор ори ) 
m [а реј и еши шш Тр өр аш ж ше е) 
ы ДАДА ДАДА А РО ү ашаа еа шы [ш ЫЕ ЕЕ Еа ае нене шш шшш [ши 
аса ЕШ EBEBBBEEBEBEBBEBESUBBEBUSEBEDSEBECEBEBEBESSESSSSSSESSSSSL | |l . 
ЕЕРЕЕ ЕЕЕ ЕГЕ ЕЕ Веща ааа | _ 
E. крш шше аши ишш и шыша 4c ------р-- 4 —- c pet 
НАТТЫ ЕЕЕ анаарар шш р к = е р те Гаји Ини 
BEE EEE EE EEE EEE EEE EEE ЕТЕ ТЕ ЕНШ ЕШ 
ШЕ ЕЕ ЕЕ Еа шшш ЕЕ жш шна ШШ И) СЕ герар а Еа а жш Т ЕТЕШ ш ваш иаа шш) | 
ЕНЕ ЕЕ На а ана 
SHES PT ere ee ee isis Sea a, 
EEE EEE e ЕН EEE 
me HH ER раса H ee Eee ee ae 
pure арыш BE [er [rm ТЕТ рит ут та 22588 CELI O omamo Jb 
e a ЖЕ 
H ыы. је eiaa ЕЕ 
= SS ee eS St eS SS SS x ‚гт ILI EL ml 
—H-44 444 8SSSSSH Se eae 
Joe а ви (ЕШ (ЕШ E] hi1 -31313 331—553 
КЕШЕГЕ ЕТЕТ Г phim SS SS SS 114—4 
pz е ps spo pne во es [ee ee Jee a | en Tre f eee т ү ет Г [= 
SR nnn Sn nS ete Sen Seen S808 Seeeeeeeue8 
OL a инини 5251 [jesse шыш 44 = 
Иши ш ши шш ЫШ к ш ыы ы аа ы ЕШ ЫЫ ГЕ BBBBZ CHEESE. SE LIT le a 
Bee tetas at st od с Емас. Е] == ЕНИН ЕЕ 
ВЕТТЕР ОЕТ СЕСЕ Е аг ЕЕЕ Еа ыа 017 
ШЕ ЕЕ Е Г) ж шшш Ба Lee Shae eee n foo [o] eor Po P || | | је] аи ај ти 
UBBBEBREBEREBEENBBEBEEBBEBBEBESBEBBNSEBEEBEBBENEBEBEZBLLÓZIZLISELINELUSSEGLZAIO- 00 
MEE m EE TEn FERESE EEE | Y 1: 
ү лл лл ү Лү Т . !|!. 
na иши ш и a шшш шаш ыша шаи шшш шш e fes fcn e ЕШ ГІГГІРРЕЕРІР шш 
па BEDE ele aaa aia) E EP [SES BBEBBEBST"SES?7) ни 
А ДЕ ДАДА ДОЗА ЕСТЕ ДАРЕ САСТАО РАВЕНИ ee 
ay mmm ME s) 3 ЕЕ ee) LS ee ae 
ШТЕТЕ ЫЕЕЕ ЕНЕ ЕЕ Е рәји 
РЕЈ ЈЕ | ја ји oa у, ү" а шш ен ка "нө терт тт ө шшш шиши 
ШЕЕСЕТ ЕГЕТЕ Е Ес ое Е еее шәрә руу и 


Figure 20. “PRs contour plot of iler elements at ГЕ = 0. 08 m fter 
detonation. 


The effect of mesh reflection, discussed in the previous section, can 
be seen in Figure 20 along the charge axis to the front of the charge nose as 
the shock wave transmits through the media at different speeds in this 
region of very poor match of element sizes. The pressure contours bend 
inward toward the charge nose in the small element region. By 0.18 ms, 


Figure 21, the mesh reflection effect in front of the charge is negligible. 
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Figure 22. Pressure contour plot of water elements at t = 0. 28 ms НЕР 
detonation. 
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Figure 23. Pressure contour plot of water elements at t = 0.38 ms after 


detonation. 
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Figure 24. Pressure contour plot of water elements at t = 0.48 ms after 
detonation. 
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Figure 25. Pressure contour plot of water elements at t = 0.58 ms after 
detonation. 


The shock front begins as a tear shape at early times emanating 
from the charge burn region and becoming nearly spherical on the charge 
axis off the nose. The above pressure contours show the nose portion of the 
evolving shock front to be nearly spherical in shape out to approximately 40 
degrees off the charge axis with a radius centered at the location of the 
undetonated charge nose. Once the burn region reaches the tail of the 
charge, the aft portion of the shock front expands to nearly spherical with a 
smaller radius of curvature than at the forward end. The above figures 
clearly illustrate the shock rapidly travelling outward ahead of the 


expanding detonation products. 
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2.  Directionality of Peak Pressure 
This subsection examines variations in the peak pressure 
developed in regions surrounding the tapered charge based upon the FEA 
conducted. Figure 19 marks the locations of the elements whose 
pressure-time histories were gathered and processed to form the plots. 


For the plots of this section, relative peak pressure P, , is defined 





as: 
Dus 
D = P. 

Where P may = maximum pressure computed for the element 


P, - maximum pressure observed for the element at the same 
standoff and nearest the charge nose axis. 
Figures accompanying the plots serve to further clarify the determination of 
P „as well as explain the geometry corresponding to each plot. 
Figure 26 corresponds to Figure 27, a plot of relative peak pressure 
as a function of 0 degrees off the charge nose axis from 8 = 0 to 90 degrees, at 


constant standoff, for three different standoffs. 


Target 


Constant R 








ee trm E Po 
Figure 26. Constant standoff, 0 degrees off charge axis. 
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Figure 27. Relative peak pressure as a function of the angle off the 
charge axis at constant standoff. 


For the three standoff distances of Figure 27, there 1s a decrease in 
pressure from the on-axis pressure as the off-axis angle increases from 0 to 
20 degrees. This decrease, ranging from 3 per cent for R = 43 inches to 5 per 
cent for R = 15 inches, may be attributed to mesh geometry and mesh 
reflections. The effect of mesh geometry was caused by the centroid of each 
element not being at the exact standoff. This effect was minimized by 


applying a first order correction to each measured pressure: 


/ 


R 
Dau Pug 
Where Pr = pressure used for plot at the nominal standoff R 


Р, = max pressure for the given element at standoff R' 


/ 
R = Jx? +y? +z’ 
x, y, z = coordinates of the element centroid. 
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After approximately 20 degrees, the curves of Figure 27 rise 
steadily with increasing angle off the charge axis. This was to be expected 
due to the tradeoff between peak pressure duration and magnitude of the 
tapered charge. As O0 increases, the target is placed nearer to the bulk of the 
mass of the charge. In a test situation gages must often be located off the 
charge axis to record the free-field pressure experienced by the target 
without being overly influenced by proximity to the target. The test designer 
often attempts to locate the gage as far off axis as possible while minimizing 
the deviation of the gage measurement from the on-axis values. Figure 27 
indicate that, for the charge and standoffs of this study, gages could be 
located up to 35 degrees off the charge axis for an error in maximum peak 
pressure measurement of less than 5 per cent. To keep the error in 
measurement of plateau duration to an acceptable level, the allowable 
off-axis angle may lie well below the level based upon peak pressure alone. 

Past 90 degrees off the front charge axis, constant standoff was 
replaced by constant distance y off the charge axis to examine the relative 
pressures experienced by elements to the side of the charge. Figure 28 
illustrates the geometry involved. P. for this case is taken at a standoff R = y 
on the axis off the charge nose as shown in the figure. Figure 29 shows 
relative peak pressure as a function of & / L off the charge axis for & / L=0 to 
] from nose to tail. As can be seen from the figure, relative peak pressure 
shows marked departure for different values of y. While the highest relative 
peak pressure occurred at €/L= 0.8, the magnitude fell from 8.4 at y = 19 


inches to 6.2 at y= 43 inches. This was likely a close-in phenomenon. Itis 
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Target 







Constant y 


Sn . Charge Axis 


Figure 28. Constant distance off the charge axis, to the side of the 
charge. 











g/l 





Figure 29. Relative peak pressure as a function of the fraction of the 
length aft of the nose at constant distance off the charge axis to the side. 


expected that, at distances further removed from the charge axis than those 
of this study, the maximum relative peak pressure to the side of the charge 
will continue to fall then reach a steady value significantly lower than 


calculated here. Of interest is the fact that the maximum relative peak 
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pressure observed for elements located to the side of the charge was found at 
E€ / L = 0.8, significantly aft of the center of gravity location at & / L = 0.6. 
This may be due to a build-up in pressure away from the center of gravity 
location in the direction of the tapered charge high explosive burn. 
Completing the examination of pressures encircling the charge 
from fore to aft, Figure 30 shows the geometry corresponding to Figure 31 
plotting the variation in relative peak pressure at a constant standoff R 
measured from the center of the tail. The angle 6 in this case is measured 
from 90 degrees off the rear charge axis toward the rear charge axis in order 
to continue proceeding in a counterclockwise direction. As in the two other 
cases, P, is measured off the charge nose, in this case at a standoff equal to 
the constant R. As in the constant R case off the charge nose, a correction 
was applied to the calculated values to account for variation in standoff of the 


element centroids. 





Target 
5 Constant R 


Charge Axis 4 R 
: a зев > —— a = 
| Po 


Figure 30. Constant R off the rear charge axis. 


Figure 31 shows a continuation of the decrease in peak pressure 


for the standoffs studied until 6 = 60 degrees, or 30 degrees off the charge tail 
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Figure 31. Relative peak pressure as a function of degrees aft of the 

charge. 
axis followed by a slight increase to the charge axis. Here again. the 
unevenness of the curves near the axis may be attributed in part to mesh 
reflections. As in the side case, the relative pressure curves, though 
exhibiting similar shape, varried significantly for different standoffs, with 
shorter standoff corresponding to higher relative peak pressure. Though of 
minor interest in nuclear simmulations, the side and rear relative pressure 
plots are of more interest to weapons and industrial high explosive designers. 
The higher pressures at the tail end of the charge agree with the use of this 
configuration for demolition and other applications. 

Of more interest in simulations is the effect of distance off axis at a 
constant standoff from the charge nose as illustrated in Figure 32. This 
configuration, somewhat easier to set up for experimental verification, 
yielded the plot of Figure 33 which provides information of similar utility as 


the constant R plot of Figure 27. 
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Figure 338. Relative peak pressure at constant distance along charge 
axis as a function of degrees off the charge axis. 


Figure 33 indicates that, for the studied charge design and stand 
off distances, the off axis distance may be varied up to 10 degrees before 
exceeding a 5 per cent error in peak pressure readings. Since the allowable 
angle to stay within a given error tolerence decreased with range, a smaller 


angle would be allowed at the longer ranges used in simulations of the 
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underwater nuclear explosion. Having concluded the FEA of the underwater 
tapered charge explosion, the next chapter outlines computer optimization of 


a simple method to determine the pressure-time history of a tapered charge. 
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V. TAPERED CHARGE DESIGN OPTIMIZATION 


This chapter outlines the coupling of a public domain computer 
optimization package ADS (Vanderplaats, 1984) and a simple tapered charge 
pressure-time generating program TAPER (Costanzo, 1991) to optimize the 
tapered charge design process. The coupling program, listed in its entirety in 
Appendix C, can be used with existing and future tapered charge 
pressure-time history generating codes to enhance the design of the tapered 
charges used to simulate underwater nuclear explosions. 

In addition to the program in Appendix C, ADS (Automated Design 
Synthesis) and a pressure-time generating subroutine are required to 


perform tapered charge design optimization. 


A. SIMPLE PRESSURE-TIME HISTORY ALGORITHM 

Simpler, less computationally intensive, computer codes than that used 
for the finite element analvsis of the previous chapter exist to predict the 
pressure profile of an underwater tapered charge explosion. These simpler 
codes are usually based upon an empirically based superposition principle. 

Because the empirically based exponential approximation of Chapter II 
works well for spherical / near spherical charges, one method of deriving the 
pressure-time history of an underwater tapered charge explosion is to 
partition the tapered charge into subsegments as shown in Figure 34. Each 
subsegment is considered to be a separate charge generating its own 


pressure-time history. (Costanzo, 1991) 
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Figure 34. Tapered charge discretized into subsegments for calculation of 
pressure-time history using summation method. 


Shocks generated by subsegments detonated after the nose subsegment 


travel in shocked media at faster speeds than preceding shocks, overtaking 


them. There is thus a summation or stacking effect of individual "wavelets" 


to determine the overall pressure-time history of the tapered charge 


underwater explosion (Costanzo. 1991): 


Where 


N А, PB: 
[= DC DOE 
i=] ^! 


P - pressure as a function of time 

N = number of waves stacked 

A; = empirical amplification factor, a function of 0; 
0; = slope angle of subsegment 

т, = R+€é; = subsegment standoff 

R = standoff from charge nose to target 

Ё = distance from nose to subsegment midpoint 
B = empirical decay constant 

t = time 

Di = charge diameter at subsegment midpoint 
66 2 subsegment length. 


Figure 35 illustrates the superposition scheme for a simple tapered charge 


discretized into large segments. 
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Figure 35. Summation of wavelets generated by charge subsegments to 
obtain resultant pressure-time history for a tapered charge. 

One algorithm, using a superposition scheme and empirical data to 
determine the pressure-time history of a tapered charge detonated 
underwater is the PASCAL program TAPER written by Fred Costanzo 
(1991). This program was converted to a FORTRAN program then to the 
FORTRAN subroutine used for tapered charge design optimization. Besides 
a basic overhaul of the input / output, a simple provision for surface cutout 


was added to the original algorithm. 


B. OPTIMIZATION PROBLEM 

The objective of the tapered charge design optimization was to 
determine charge geometry, as depicted in Figure 8, and standoff required to 
develop a pressure-time history most nearly matching a desired pressure- 
time history. The average square root of the sum of the squares of the 


differences between the computed and desired pressures at each time 
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increment was selected to quantify the difference between the optimal design 
and desired pressure profiles. Constraints restricted the charge to increasing 
diameters and placed upper and lower limits on charge weight, charge 
dimensions, and standoff. Only designs within these limits were allowed. 


For a tapered charge of three segments, the optimization problem became: 





minimize: |2 K к= N 
subject to: lin <L < 20 ft, 1= 1 0 3 
0.7/5 1n « d, « 10 1n, 1= 1104 
5 ft « R « 30 ft, 
25 lbm < W < 125 Ibm. 
= а еа, 
where N = number of computed pressure-time history points 


P, = desired pressure at a particular time 
p, = computed pressure at a particular time 
L = charge segment length 
d, = charge joint diameter 
R = standoff 
W = charge weight. 
The main FORTRAN program DTAPOPT was written to accomplish the 
tapered charge design optimization using, as described previously, the ADS 


optimization package and the subroutine based upon the TAPER program. 


C. OPTIMIZATION PROGRAM 

In addition to DTAPOPT, three subroutines for use in conjunction with 
the main program were also written. The first, DCOMPAR, computes the 
square root of the average sum of the squares of the pressure differences. 


The second, DPTGEN, interpolates to find the desired pressure 
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corresponding to the design pressure at a given time. This subroutine 
enables the comparison of the desired and design pressures at the same 
times. The third subroutine, DTAPWT, computes the weight of the tapered 
charge for evaluation of the weight constraint. 
1. Required Program Input 

DTAPOPT requires input from two files and the kevboard. Figure 
36 shows a sample of a tapin.dat file used to input initial charge geometry, 
standoff, nominal length of pressure-time history to be computed, and surface 
cutoff time. Decimal alignment and horizontal placement on the lines is 
optional with one number per line only. Number translations have been 


written into the figure. 


number of charge segments 
length of first segment (ft) 
length of second segment (ft) 
length of third segment (ft) 

first joint, or nose, diameter (in) 
second joint diameter (in) 

third joint diameter (in) 

fourth joint, tail, diameter (їп) 
standoff (U 

1.12 nominal length of computed time history (ms) 
1.095 surface cutoff time (ms) 


3 

Ge 
ШІ 
S 
l. 
D 
24 
2 
l 





Figure 36. Sample tapin.dat file. Italics, not part of the file, indicate 
the meaning of each number. 

Figure 37 shows a sample of the second input file, profin.dat, used 
to input the desired pressure profile for program DTAPOPT. Up to 1001 data 
points may be entered without program modification. The first time entered 
must be at time zero, and the last time must be greater than the nominal 


length of computed time history entered in tapin.dat. Decimal alignment 
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and horizontal placement on the lines is optional. Only one number can be 
used on the first line, two on the rest. Number translations added to the 


figure are shown in italics. 





ndes, total number of desired pressure-time data points minus 1 





0.0 0.0  tdes(0) in ms--must be zero. pdes(0) in psi 
0.04 1950.0  tdes(l) pdes(1) 
0.119 1755.0  tdes(2 pdes(2) 
0.278 1560.0  tdes(3) pdes(3) 
0.417 1365.0  Ides( 1) pdes( 1) 
8855501370.0 . fdes(3) pdes(5) 
0.834 975.0 tdes(6) pdes(6) 
1.072 780.0  tdes(7) pdes(7) 

| ШО 120 0.0 14е<(%) раез(8) 

| 3 0 0.0  tdes(9) pdes(9) 





tdes(ndes) must be greater than nominal length of p-t history entered in tapin.dat 





Figure 37. Portion of a sample profin.dat file. Italics, not part of the 
| file, have been added to indicate what each number represents. 


Three integers comprise required keyboard input. These numbers 
control the optimization method used by ADS. Vanderplaats (1985) provides 
more detailed instruction on method selection. Vanderplaats (1984) provides 
the theory behind the methods. The three numbers consist of any one 
number from each of three groups. The first number may be any of the 


following to determining the optimization strategy: 


First 
Number Optimization Strategy 
0 Go directly to the optimizer 
6 Sequential Linear Programming 
7 Method of Centers (Design must be feasible) 
8 Sequential Quadratic Programming 
9 Sequential Convex Programming. 


49 


The second required input number may be either of the following to 


determine the optimizer: 


Second 
Number Optimizer 
4 Method of Feasible Directions 
5 Modified Method of Feasible Directions. 


The last input number. one of the following. selects the one-dimensional 


search option to be used: 


Third 
Number One-Dimensional Search 
2 Golden Section Method 
6 Golden Section Method Plus 
Polynomial Interpolation 
T Bounded Polynomial Interpolation 
8 Unbounded Polvnomial Interpolation. 


2. Program Output 

Output from DTAPOPT includes one screen summarizing the 
optimization and an output file containing the design and interpolated 
desired pressure-time histories. 

Figure 38 shows a sample screen from a run of DTAPOPT on a 
personal computer. The figure shows the execution commands followed by 
prompts for the three inputs, in this case the user selected the combination 
8-5-7 for the strategy, optimizer, and search. The output summary then lists 
the initial and final charge designs as well as the average square root of the 
pressure differences for each case and the number of calls to the pressure- 


time generating routine. 
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C:\P&C\FOR\OPT: dtapopt 
ENTER ISTRAT.IOPT.IONED: 8 5 7 


INITIAL DESIGN 
LENGTHS = 1.000 1.000 5.000 ЕТ 
DIAMETERS = 1.00 2.00 4.00 5.00 IN 
CHARGE WEIGHT = 66.0 LB 
RANGE = ОФСЕТ 


SQRT OF AVG OF (Pdesign-Pdesired)2 - 424.3 PSI 


FINAL DESIGN 
LENGTHS = 0.083 0.083 4.825 FT 
DIAMETERS = 2-50 SS 4.28 4.52 IN 
CHARGE WEIGHT = 55.7 LB 
RANGE = ATI 
SQRT OF AVG OF (Pdesign-Pdesired) 2 = 240.3 PSI 


CALLS TO P-T HISTORY GENERATOR - 158 





C:\P&C\FOR\OPT: 


Figure 38. Sample screen output from DTAPOPT. 


The output screen of Figure 38 represents an early step in the 
design process. The next would be to use this final design to input a more 
refined initial design, then run the optimization program again. For test 
optimizations, the square root of the average pressure difference squared was 
well below 100 for the "best" final design. 

Figure 39 shows a portion of a sample DTAPOPT output file. The 
first two columns contain the calculated time and pressure, the third column 
contains the interpolated values of the desired pressure profile input. The 
data are in free format to retain maximum precision, sacrificing the 
readability of formatting. This output file may be readily used to create plots 
comparing the initial and final designs as was done in the next section of this 


chapter. 
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0.00000000000000000E-01 0.00000000000000000E-01 0.00000000000000000E-01 
2.81022269002981981bE-02  1639.32749805998742000 1405.11134501490983000 
5.62044538005963962Е-02 1675.89437624031484000 1987.35262142392480000 


0.47773785730506935 1628.10858635031309000 1658.35094063994575000 
0.50584008420536752 1614.83610147863396000 1636.41749525434739000 
0.53394231110566581 1602.96393658242482000 1614.4840498687485 7000 


1.06788462221133162 900.78633515776994000 899.17511224684437800 
1.09598684911162980 0.00000000000000000E-01 235.12876911529235700 





Figure 39. Sample excerpted from an output data file generated by 
DTAPOPT. 


D. OPTIMIZATION RESULTS 

The design space using the subroutine adapted from program TAPER 
proved to be fraught with local minima attributable discontinuities resulting 
from integer changes in the number of waves stacked. Several runs of 
DTAPOPT using different optimization methods for a partcular initial 
design resulted in an improved design unless the initial design was optimal. 
The improved design was then used as the initial design for further 
optimization. This process was repeated from four to seven times until 
further optimizations failed to improve the design an appreciable amount. 
Each run of DTAPOPT took approximately two minutes on a personal 
computer of modest, 386SX, capacity. Total time to perform a tapered charge 
design optimization was from one to two hours. 

Of the desired pressure profiles and initial designs tested, input 
parameter combinations 0-5-7, 8-5-7, and 9-5-7, usually produced the most 


improved design with the fewest calls to the pressure-time history generating 
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subroutine. These input numbers translate, as shown in Section V.C.1 to 
direct optimization, sequential quadratic programming, or sequential convex 
programming strategy combined with the modified method of feasible 
directions optimizer and а bounded _ polynomial interpolation 
one-dimensional search. For the best designs found, as mentioned in Section 
V.C.2, the square root of the average difference between desired and design 
pressures was well below 100 psi. Figure 40 shows a comparison of the 
pressure profile generated by an optimized design from program DTAPOPT 


and the corresponding desired pressure profile. 

















DESIRED ———— 
DESIGN —1— 


| 
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Figure 40. Pressure profile resulting from optimization compared with 
the corresponding desired pressure profile. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The boundary element / finite element analysis of the side on attack of a 
simple cvlinder by nuclear and conventional type pressure profiles resulted 
in substantially different structural responses. It is therefore necessary to 
model the duration of the pressure plateau to adequately simulate structural 
response to an underwater nuclear explosion. 

The finite element analysis of an underwater tapered charge explosion 
provided insight to the early time propagation of the shock wave and the 
directional variations in peak pressures developed in the surrounding media. 
There is need for future research, including sensitivity studies of the 
artificial viscosity coefficient and the time integration step, to obtain less 
oscillatory results to provide tapered charge pressure plateau duration 
information. Additionally. the computationally intensive method used lends 
itself more readily to supercomputer use where the mesh size may be 
extended far enough away from the charge form comparison with test data. 

Coupling of an optimization routine with a tapered charge pressure 
profile generating routine provides a tool which can be used to more 
efficiently design tapered charges used to simulate underwater nuclear 
explosions. The program written for this study may be used with existing 


and future pressure-time history codes. 
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APPENDIX A: TAPERED CHARGE FEM INPUT FILE 


Following 1s a complete listing of the input file used to perform FEM 


analvsis on the tapered charge and water model of Chapter IV: 


GOGO2 


cu up for DYNA3D, 
10 


integrate to time 2000, TAURUS data dump interval 


С high speed printer dump interval 29999. 


dn3d vec term 2000 


plti 10.0 prti 29999.0 


c Scale x, y, z axes from inches to cm. 


ШЕЯ 2 54 у<са 2.54 zsca 2.54 


с Define material 1, 
оте топа 1топ vel O. 


type 8--High explosive burn. HBX-1: 
731 cm/us, Chapman-Jouget pres 0.2204 Mbar. 


EE eusrty 1./12 gm/cc. 


пас 1 type 8 d .731 рсј .2204 ro 1.712 


c Equation of state 


for charge: 2--JWL: 


КЕ 2007. 5 .08513, х1 4.90. r2 1.10, omega .20, eO .103 Mbar-cc/cc. 


E03 7.5807 b .08513 r1 4.90 r2 1.10 omega 0.20 eO 0.105 endmat 


c Define material 2, 


РИЕЕЦПСОЦС pressure . 


type 9--Null Material. Water: 
1l psi (6.89e-9 Mbar), density @ 4 deg C 1.000 gm/cc. 


mat 2 type 9 pc 6.89e-9 ro 1.000 


c Equation of state 
с sound speed .142 
c Gruneisen gamma . 


eos 4 sp .142 sl 2. 


for water: 4--Gruneisen: 
спиш О degs C. sl 2.56. s2 -1.986. s3 .2268, 
4934, first order vol cor'n to gamma 1.3937. 


56 s2 -1.986 s3 0.2268 gamma .4934 sa 1.3937 endmat 


c Define two symmetry planes by defining a point and a normal vector 


c for each plane. 


c included in that 


plane 2 
Коо 100 001 
ЕО о 010 .001 


с Define detonation 


@etp 1 point О О О; 


Any point within .001 cm of a symmetry plane will be 
plane’s definition. 

symm 

symm 


Points in HBX.l. 
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c Define rear and front water cylinders. cone and plane surfaces for 
c charge nose (part 1) and transition water (part 2). 


sd 1*en2p ОБОИ ТИТ 05525 ОРОЛИ 5125 -4.0 
sd 5 plan 0 OUEST 
pow 


58 4 plan 0O O -4 O 


c Part 1: small charge cone. 


2а UIS СИМЕ 

sti '«Lus35 1.3 d 
sfi l1 -3wel-32 20989992 
dod — 32297 

mate 1l 
end 


c Part 2: water transition from small charge cone to залася 


4 
4 
4 i EN. 

Str--L 517-5 ЭКСПЕ sd 3 
E -l-5 2 2 ва 
3 
5 


c Define cone and plane surfs for med charge and transition (pts 3 and 


sd 8 plan 0 0 - 


бада 
СОС спри ОИЕ eo ae -4.0 2.0625 -8.0 


c Part 3: medium charge cone. 
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ПО 3: 103: : 

BENI -5. -1-3: : sa 9 
2E! -3: -1 -3: -1 -1: sd 4 
СИ Э -1 -3; -2 -2; sd 8 
NN 2232 

UM] 322 

mate 1l 
end 


СШ ОШ Л medium water transition. 


Smart | 3 5 7/9: 





Off ~ 
гг 


eee Ob Aa cs vr 


c Define cones and plane for large charge and transition (pts 5 and 6). 


за 11 спрроо 0 001 2.0625 59-0 E6875 -60.0 
sd 13 plan 0 0 -60 001 


c Part 5: large charge cone. 


БН -! -3: -1 -3; ; sd 11 

sfi |І -3; -1 -3: -1 -1; sd 8 
EN 1] -5: -1 -3; -2 -2; sd 13 
ШЕГІ 23) 

НА 322 

mate 1l 

end 


c Part 6: large water transition. 


ОСО А Аи 
sfi =7-=4% -2 ЕСЕ 
sfi. -L.-5 : -1755 Tem с 


sfi -І -5 :-І1 57 -2--2777 СЕБЕ 
d РОЗА 

алта” о 

mate 2 
end 


c Define cylinders for front and rear water cylinders and transitions. 
sd ссу Оо Оаа) 200575 
c Part /: rear water cylinder. 


start 
Э. 
5 


* 
, 


NW WR 


4: 

10 1 

-1 0 1 

-60.0 -106.0 

ато з Јово 

БЕ) (1 ЭЕ ЕЭ сахех 

sti -1-37 51 ау = = sd ES 


IP eect 


dq I 2$ 
db DEUS 397 2 
nro ШІ 009272 0802 
mate 2 

end 


c Part 8: rear water transition. 


0 
4 
5. eL 55 еј = А ЛЕО 
5 
3 


(іле ӘБЕУ 009091922 
mate 2 
end 


c Define cones and plane for water cone and transition (parts 9 and 10). 


sd 5 cn2p 00 0 OO OF 5625 0.0 DESEAS A 
sd 7 plan 00 4 82207! 


C Part 9: front water conc. 


start 


tar 
LONS 
10957 


1 
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I5. -I1 E. ' sd 5 
Peet -3; -1 -3; -1 -1; sd 3 
sfi -l -3; -1 -3; -2 -2:; sd 7 


DENEN 2 -4; -2 -4: ; sd 5 

sfi -1l -5; -1 -5; -1 -1:; sd 3 
EE] -5: -1 -5; -2 -2: sd 7 
БИ !І 352 

КИН 532 

mate 2 
end 





Себе cylinder for front water. parts 11 апа 12. 
БИЕ СІП 000 001 1875129 


c Part 11: front water cylinder. 
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(4 
а 1 2 08595: 19200920 9 


аб ЈА | ae? Sede 

ЕГ 1 -57 Бр е Ва сае 
ce PNIS 
ІР ІІ Ө 
пр ОО 2 ПО? 
mate 2 

end 


Sle 
5 2 
3087 


c Part 13: top water. 


start 

] 25 

Bo 

LEA SA: 
0.0 46.0 
Оов 
-106.0 -60.0 
ЙҮ Б 020g» 
nrb 02250050 
nib 00210) 
перо О В 
mate 2 

end 


c Part 14: side water. 


start 
22. 
es: 
TED 54-277: 
40 46.0 
О Он у 
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APPENDIX B: GRUNEISEN GAMMA APPROXIMATION 


The seventh order approximation for the Gruneisen gamma as a 
function of the specific volume v developed by Gurtman, Kirsch, and 


Hastings (1971) is: 


уү(у) = ао+ају +аоу“ +. · ·+ату' 


where ao = 2,366.6324 

a, = -22,669.420 

ao = 91, 259.368 

аз = – 200, 175.85 

20— 255-585. TI 

a5 — -196, 872.84 

ас = 81, 850.023 

ат = 14. 842.530 

у = Gruneisen gamma, dimensionless 
cm. 


em 





v = specific volume in 


By definition. 


1 


Vo 
ЕНІ a. 
у папа C 


= = 
util w+il 





giving 


Substituting l/(jt* 1) for v into the Gurtman equation for u = 0 to 0.8 in 
.001 increments, then using a least squares linear fit with the same intercept 


resulted in the following linear equation for y(t): 


y(u) = yo + ay = 0.4934 + 1.3937 
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Figure 41 is a plot of y(u) from the Gurtman equation and the lnear 


approximation. 


Gruncisen Gamma 











9=.4934+1.393/т 





gamma from 
Gurtman 








== 


06 07 





Figure 41. Comparison of y(u) from Gurtman, Kirsch, and Hastings 
equation with the linearized equation for yu). 
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APPENDIX C: FORTRAN PROGRAM 


This appendix contains a complete listing of the FORTRAN main 
program DTAPOPT and its supporting subroutines DCOMPAR, DPTGEN, 
and DTAPWT, written to optimize tapered charge design. A separate 
subroutine to generate pressure-time histories given tapered charge 
geometry and standoff distance is required, as well as the optimization 
package ADS. 

PROGRAMMING NOTE: The program DTAPOPT and associated 
subroutines were written in double precision FORTRAN. The public domain 
version of ADS is written in single precision. ADS was converted to double 
precision for use with DTAPOPT. DTAPOPT and its three subroutines can 
be simply converted to single precision by removing the IMPLICIT NONE 
and DOUBLE PRECISION statements from the codes. 


I. PROGRAM DTAPOPT 


DTAPOPT 

This FORTRAN program combined with subroutines DCOMPAR, 
DPTGEN, and DTAPWT, is designed to optimize tapered charge design 
using the public domain optimization package ADS (converted to 
DOUBLE PRECISION by the author) 

--When coupled with a separate subroutine (not included) 
to calculate the pressure-time history of a three-segment tapered 
charge. 

EM to ten charge segments may be used with appropriate 
modifications to the input files. 


PRECISION: DOUBLE 
MEUT FILES: TAPIN.DAT Initial Design. 
PROFIN.DAT Desired P-T History. 
INTERRACTIVE INPUT: Optimization options. 
OUTPUT FILE: Named in TAPIN.DAT, plot data for computed 
time vs computed pressure and interpolated pressure. 
SCREEN OUTPUT: Starting and Optimized Designs. 
REQUIRED SUBROUTINES: DTAPWT Computes charge weight. 
DPTGEN Linearly interpolates to give desired 


(0 Q (0 0 оосо ооо О ОООО ОШОК ОШО 
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pressure at each time output by P-T history generator. 

DCOMPAR Computes the square root of the 
average of the square of the difference between the computed and 
interpolated desired pressure at each time. 

ADS Package of subroutines which perform 
the optimization. 

P-T HISTORY GENERATOR. 


AUTHOR: William Earl Miller II 
MOST RECENT UPDATE: 6/16/92 
REFERENCES: (1) Vanderplaats, G. N.. ADS - A FORTRAN 


Program for Automated Design Synthesis, Version 1.10. program 
instructions, Naval Postgraduate School, Montery, California, 
May. 1985. 
(2) Vanderplaats, G. N., Numerical 
Optimization Techniques for Engineering Design: With Applica- 
tions, McGraw-Hill Publishing Company, 1984. 
TEST P-T GENERATOR: The subroutine SDTAPER., was converted to 
FORTRAN by the author, based upon the PASCAL program TAPER version 
7/27/89 by b. А Gostanzo, 
DATA INPUT FILES--Two required: 
TAPIN.DAT contains 12 lines. one value per line which are read 
into the program. The integer and decimal data need only be in 
a format suitable for list-directed input assignment to INTEGER 
and DOUBLE PRECISION data types respectively. The character data 
must be in proper form for a DOS file name. 


LINE READ TO DATA 
NUMBER DESCRIPTION VARIABLE ТҮРЕ 
1 number of charge segments NTAP integer 
2 length of charge segment (ft) TAP(1) decimal 
D Қ TAP(2) decimal 
4 E ТАРС) decimal 
5 diameter of charge segment (in) DIAM(1) decimal 
6 DIAM(2) decimal 
7 ч DIAM(3) decimal 
8 M X DIAM(A) decimal 
9 standoff (IE) RANGE decimal 
10 nominal length of time history (ms) TLEN decimal 
used by P-T generating subroutine 
ІШІ pressure cutout time (ms) TCO decimal 
used by P-T generating subroutine 
12 name of output data file to be NAMFIL character 
created 


PROFIN.DAT contains NDES+2 lines, one integer value (NDES) on the 
first line, two decimal values on each of the remaining lines. 
Times must be in ascending order starting with 0.0 and extending 
to a time greater than TLEN entered in TAPIN.DAT above. 


LINE READ TO DATA 
NUMBER DESCRIRETON VARIABLE TYPE 
l total number of desired pressure- NDES integer 
time history data pairs minus one 
2 first time first pressure TDESCO) PDESCUD decimal 
(must be zero) 


3 second time second pressure TDESCIO PDESQORB) decimal 


NDES42 NDES time NDES pressure TDES(NDES) PDES(NDES) decimal 


INTERRACTIVE KEYBOARD INPUT--three integers required 
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These integers control the optimization method used by ADS. 

See Ref 1 for more complete instructions, Ref 2 for theory. 
Combinations O 5 7, 8 5 7, and 9 5 7 are recommended. Others may 
work well in a given design space. Using the initial design. three 
or four runs with different combinations should result in a most 
improved design which can then be used as the initial design for 
further optimization. The three integers are: 


ENTRY READ TO DATA 
ORDER DESCRIPTION VARIABLE TYPE 
first optimization strategy used by ADS ISTRAT integer 


O Go directly to the optimizer 
6 Sequential Linear Programming 
7 Method of Centers (Design must be feasible) 
8 Sequential Quadratic Programming 
9 Sequential Convex Programming 
second optimizer to be used by ADS LOPE integer 
4 Method of Feasible Directions 
5 Modified Method of Feasible Directions 
third one dimensional search options IONED integer 
5 Golden Section Method 
6 Golden Section Method Plus Polynomial Interpolation 
7 Bounded Polynomial Interpolation 
8 Unbounded Polynomial Interpolation 
OUTPUT FILE: Named in TAPIN.DAT. consists of NTIME+1 rows of 
three columns. Format is list-directed from DP variables. 


COLUMN WRITTEN FROM 
NUMBER DESCRIPTION VARIABLE 
one times output from P-T GENERATOR ТІМЕ(І), I-O,NTIME 
two P's from P-T GEN'R at each time PRESSCI), I-O,NTIME 


three int. P’s for each time from DPTGEN PCOMPAR(I1), I=0,NTIME 
SCREEN OUTPUT: Outputs initial and final design values plut the 
number of calls to the P-T History Generator. 

INITIAL AND FINAL VALUES 


DESCRIPTION VARIABLE 
lengths of tapered charge segments (ft) TAP(I),I-1.NTAP 
diameters of tapered charge (in) РТАМ(1Т), 120. МТАР 
charge weight (lb) WEIGHT 
standoff (Тү) RANGE 
sqrt of average sq diff of Pdesign-Pdesired OBJ 
OPTIMIZATION EFFICIENCY 
number of calls to P-T generator NCALLS 
VARIABLES 
NAME DATA | COMMON ASSIGNED USED DESCRIPTION ADS 
ТҮРЕ BY | BY ae 
BEND), DP - ADS ADS only Constr grads Y 
I-1,NRA;,J-1,NCOLA 
БЕ( 21) DE - ADS ADS only Obj grads Y 
DIAM(I), DP OPTDPA TAPIN.DAT,MAIN MAIN, Charge diams N 
I=] ,NTAP+1 PTGEN’R 
BI), DP - MAIN ADS Constraints ү 
I=1 , NCON 
| 1 - MAIN MAIN Local indexing N 
ПОСТО, i - ADS ADS only Gradient ID ү 
I-1,NCON 
IDG(I), I - MAIN ADS Constraint ID ү 
1=] ,МСОМ 
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930200000390) 70 OQ QOO (1) 000 (Q0 0 0 0D O0 0 O0. 0 0 СУ СУ Qe Gee? OOO} OO 0} O20 Oo) Ко ОУН СУ (220 Оо Оо СеО ОВИЕ О 


IGRAD 
INFO 

IONED 

IOPT 
IPRINT 
ISTRAT 
IWK(I), 
I-1,NRIWK 

NAMFIL 

NCALLS 

NCOLA 

NCON 

NDES 

NDV 

МСТ 

МВА 

NRIWK 

NRWK 

NTAP 

NTIME 

OBJ 

PCOMPAR(1), 
I-0,NDES 

PDES(I) 
I-0,NDES 

PRESS(I) 
I-0,NTIME 

RANGE 

SUMSQ 

TAP(I) 
I-1,NTAP 

TCO 

IDEST 
I-0.NDES 

TIME(I) 
I-0,NTIME 

TLEN 

VLB(I), 
I=1 , NDV 

VUB (IDE 
I-1,NDV 

WEIGHT 

WK(I), 
I-1.NRWK 
І-1,Мру 


SUBROUTINES 


Е IM eS I ЕЗІ ЕЗІ 


C 


I— I I HA нч нч нч нч нч EA A 5 


DP 


FORDPTGEN 


OPTI 
OPTI 


FORDPTGEN 
FORDPTGEN 
OPTDPA 


OPTDP 


FORDCOMPAR 


OPTDEA 


OPTDP 
FORDPTGEN 


OPTDPA 


OPTDE 


PASS 


MAIN ADS Grad Calc Ctrl 
MAIN,ADS MAIN,ADS Prog Flow Ctrl 
user ADS One-D Search 
MAIN ADS Optimizer 

MAIN ADS ADS Print Opts 
user ADS Strategy 
ADS .MAIN ADS Work Array 
TAPIN.DAT MAIN Output File 
MAIN MAIN PTGEN'RWecalls 
MAIN ADS No. A columns 
MAIN ADS No. constraints 
PROFIN.DAT DPTGEN Des'd plot sts-l 
MAIN ADS No. Design vars 
ADS ADS Gradient ctrl 
MAIN ADS A rows 
MAIN ADS IWK Dimension 
MAIN ADS WK Dimension 


TAPIN.DAT PTGEN'R,MAIN No. Chg Segs 
PTGEN'R DPTGEN,DCOMPAR Plt sts-l 


MAIN MAIN ADS Obj Func Val 
DPTGEN DCOMPAR Interp’d Press 
PROFIN. DAT DPTGEN Desired Press 
PTGEN'R PCOMPAR,MAIN Calc Press 
TAPIN.DAT,MAIN PTGEN'R Standoff 
DCOMPAR MAIN Press var'n 


TAPIN.DAT,MAIN PTGEN'R,MAIN Chg Seg Lth 


Surf CO Time 
Des'd PT Time 


PTGEN'R 
DPTGEN 


TAPIN.DAT 
PROFIN.DAT 


PTGEN'R DPTGEN,DCOMPAR,MAIN Calc time 


TAPIN.DAT PTGEN'R Lngth PT hist 
MAIN ADS L lim on DV 
MAIN ADS U lim on DV 
DTAPWT MAIN Charge Weight 

ADS, MAIN ADS Work Array 
MAIN, ADS ADS Des Vas 


ADS(INFO,ISTRAT,IOPT,IONED,IPRINT,IGRAD,NDV,NCON.X,VLB, VUB,OBJ, 


G,IDG,NGT,IC,DF,A,NRA,NCOLA,WK , NRWK , IWK, NRIWK ) 


See Ref 1 for more complete instructions, Ref 2 for theory. 


Simply, ADS inputs design variables. constraints, and the objective 


function for an initial design, then modifies that design, 
requesting the corresponding objective and constraint values from 
the calling program. 
ARGUMENT VARIABLES 


A(NRA,NCOLA) DP Array of constraint grads. 
DF (NDV+1 ) 


G(NCON) 


ADS use only here. 


DP Array of objective gradients. ADS use only here. 
DP Array of constraints for current design in X 
G(1&2): 25#<WEIGHT<125#; G(3-5): DIAM(1)<DIAM(2)<DIAM(3)<DIAM(4) 
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IC(NGT) I Array of constraint ID's. NA this program 
IDG(NCON) I Array ID'ing type of constraints: O for nonlin ineg 
IGRAD I =0 for ADS calculate gradients using FD. 

INFO I ADS flow control parameter. 

IONED I5 6.745. See input section. 

ШОРТ I =4,5: See input section. 

IPRINT І =0000 FOR NO ADS PRINTOUT 

ISTRAT I =0.6,7,8,9: See input section 

IWK(NRIWK) I Stores ADS I vars. Some modifiable. 

NCOLA І Dimensioned columns of A, min NDV+1. 

NCON I Number of constraints in C. 

NDV I Number of design variables in X. 

NGT I Returned to by ADS for gradients. O this program. 
NRA I Dimensioned A rows: at least NDV+1. 

NRIWK I Est: 200+NDV+NCON+N+MAX(N, 2*NDV) .N=MAX(NDV,NCOLA) 
NRWK I Est: 500+10*(NDV+NCON )+NCOLA* (NCOLA+3 )+N*(N/2 )+1 
OBJ DP Objective function value, SUMSQ provided by DCOMPAR 
VLB(NDV+1 ) DP Array of design variable lower bounds, indices as X 
VUB(NDV+1 ) DP Array of design variable upper bounds. indices as X 


WK DP Array for ADS double precision variables. 


X(NDV+1 ) DP Array of design vars. Assigned by input, then ADS. 
X(I)=TAP(I), I=1,3; X(4+1)=DIAM(I). I=0,3 >; X(8)=RANGE 
ERPCI)-1" VLB(4+I1 )=0.75" VEBCS )=5’ 
VUB(I)=10' МОВ(4+1)=10" VUB(8)=30' 

DCOMPAR 

Communicates via COMMONS: OPTI, OPTDPA, FORDPTGEN, FORDCOMPAR 
Input variables: NTIME: PCOMPAR(I), PRESS(I), I=0.NTIME 
Output variable:  SUMSQ 


DPTGEN 
Communicates via COMMONS: OPTI, 
Input variables: NTIME; TIME(I), I=0.NTIME 
NDES РОБУ) РОБУ), 

Output variables: PCOMPAR(I). I=1.NTIME 
DTAPWT 

Communicates via COMMONS: OPTI, OPTDPA, PASS 

Input variables: ҚТАР: DIAM(I),I=0O,NTAP, TAP(I),I=1,NTAP 
Output variables: МІЕСНТ 
PRESSURE-TIME HISTORY GENERATOR 

Communicates via COMMONS: OPTI, OPTDP, OPTDPA 
Input variables: NTAP; RANGE; TLEN; TCO; TAP(I), 


1=0,МТАР. 
pU variables:  NTIME; TIME(I), PRESS(I), I-O,NTIME 


? 
ххх» УКУК КККК АКА КЛАСКА КУКУК УСК ЖСК АСУ УАУ АА 


OPTDPA, FORDPTGEN 


I=0 . NDES 


I=1,NTAP; DIAM(I), 


PROGRAM DTAPOPT 
IMPLICIT NONE 

SPECIFICATIONS FOR SUB VARIABLES 
COMMON/OPTI/NTAP ,NTIME 
INTEGER МТАР,МТІМЕ 
COMMON /OPTDP/ RANGE , TLEN , TCO 
DOUBLE PRECISION RANGE,TLEN,TCO 
COMMON/OPTDPA/ TAP, DIAM, TIME, PRESS 
DOUBLE PRECISION TAP(10),DIAM(0:10),TIME(0:1000) , PRESS(0:1000) 
COMMON /PASS / WEIGHT 
DOUBLE PRECISION WEIGHT 
COMMON /FORDPTGEN/NDES , TDES, 
INTEGER NDES 
DOUBLE PRECISION TDES(0:1000) , PDES(0:1000) , PCOMPAR(0:1000) 
COMMON /FORDCOMPAR/SUMSQ 
DOUBLE PRECISION  SUMSQ 


ЕВЕ; РСОМРАК 


AMAA 
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SPECIFICATIONS FOR ADS VARIABLES 

up to 20 X's, 100 С'в, ЗО соп ЕЕЕ 
INTEGER IWK(2000),IDC(100),1C(30) 
DOUBLE PRECISION X(21),VLB(21),VUB(21),G(100),DF(21), A(21,30), 


c WK(10000) 


INTEGER NRA,NCOLA,NRWK,NRIWK, IGRAD.NDV,NCON, ISTRAT, IPRINT.IOPT, 


c IONED, INFO, NGT 


DOUBLE PRECISION OBJ 
SPECIFICATIONS FOR SUBROUTINES 
EXTERNAL SDTAPER 
EXTERNAL DTAPWT 
EXTERNAL DPTGEN 
EXTERNAL DCOMPAR 
EXTERNAL ADS 
SPECIFICATIONS FOR LOCAL VARIABLES 
INTEGER I,NCALLS 
CHARACTER*12 NAMFIL 
BEGIN EXECUTION 
NCALLS=0 
ADS ARRAY DIMENSIONS 
NRA=21 
NCOLA=30 
NRWK=10000 
NRIWK=2000 
ADS PARAMETERS 
(no user provided gradients, 5 constraints) 
(f design variables determined in initial design section) 
IGRAD-O 
NCON=5 
INPUT INITIAL DESIGN, OUTPUT FILE, 
BOUND ON DESIGN VARIABLES 
OPEN(88,FILE-'TAPIN.DAT', STATUS-'OLD') 
READ (88 ,*)NTAP 
NDV-2* (NTAP*1) 
1" « length of segment « 20' 
DOS99 I-I.NTAP . 
READ( 88 ,*)TAP( I) 
XCI)-TAP(I) 
VLB(I)-1.0D0/12.0D0 
VUB(I)-20.0D0 
CONTINUE 
0,75" « joint diameter <= PO 
DO 100 I=1,NTAP+1 
READ( 88 ,*)DIAM(I-1) 
Х(МТАР+І )=рІАМ(І-1) 
УІВ(МТАР+І)=0.75ро 
VUB(NTAP+I )=10.0DO 
CONTINUE 
5’ < standoff <50/' 
READ (88 ,* )RANGE 
X(2*NTAP+2 )=RANGE 
VLB( 2*NTAP+2 )=15 .ODO 
VUB( 2*NTAP+2 )=30.0D0 
INPUT LENGTH OF TIME RECORD PICs 
READ (88 , *) TLEN 
READ(88,*)TCO 
INPUT DESIRED PRESSURE ER Gr Men 
READ (88, 1000)NAMFIL 
CLOSE(88) 
OPEN( 89, FILE=NAMFIL) 
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1000 FORMAT(A12 ) 
EDHNCSS, FILE-'PROFIN.DAT',STATUS-'OLD') 
READ(88 ,*)NDES 
DO 101 I=0,NDES 
BESDCS8, *)TDES(I).PDES(I) 
101 CONTINUE 
CLOSE(88) 
C ID FIVE NONLINEAR CONSTRAINTS 
IDG(1)=0 
IDG(2)=0 
IDG(3)=0 
IDG(4 )=0 
IDG(5)=0 
С INPUT 
C (no ADS print, interractive,optimizer and search options) 
IPRINT=0000 
BEMNIS.' ENTER ISTRAT,IOPT,IONED: ' 
EEAD*.ISTRAT,IOPT,IONED 
C OPTIMIZE 
C (0 no override, -2 default override) 
INFO--2 
ШО ADS(INFO,ISTRAT,IOPT,IONED,IPRINT.IGRAD,NDV,NCON,X.VLB, 
C VUB,OBJ,G, IDG, NGOT. IC, DF, A. NRA, NCOLA, WK, NRWK, IWK , NRIWK) 
IFS TO CONTROL FLOW USING ADS INFO 





IF(INFO.EQ.1) THEN 
EVALUATE OBJECTIVE AND CONSTRAINTS 
ADS VARIABLES TO SDTAPER INPUT 


Су С? 


ро 102 І-1.МТАР 
TAP(I)eX(I) 
102 CONTINUE 
DO 103 I=1,NTAP+1 
DIAM(I-1)=X(NTAP+I) 
103 CONTINUE 
RANGE=X( 2*NTAP+2 ) 
С CALCULATE PRESSURE PROFILE 
CALL SDTAPER 
NCALLS=NCALLS+1 


C CALCULATE WEIGHT 
CALL DTAPWT 

C CORRESPONDING DESIRED PRESSURES 
CALL DPTGEN 

с EVALUATE OBJECTIVE FUNCTION 
CALL DCOMPAR 
OBJ=SUMSQ 

C PRINT INITIAL DESIGN DATA TO SCREEN 
IF(NCALLS.EQ.1) THEN 

PRINT*, ' INITIAL DESIGN' 


PRINT1001,(TAP(I),I-1,NTAP) 
PRINT1002, (DIAM(I),I-0,NTAP) 
PRINT1003,WEIGHT 
PRINT1005,RANGE 
PRINT1004 , OBJ 
ENDIF 
C EVALUATE CONSTRAINTS (254 « W « 1254) 
G(1)2(25.0D0-WEIGHT) 
С(2)=(МЕІСНТ-75.0р0) 
C EVALUATE CONSTRAINTS (D1<D2<D3<D4) 
G3) =X(4)=X(5) 
G(4)=X(5)-X(6) 
Бр Е) ЬУ) 
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GO TO 10 


ELSE 
IFC INFO.EQ.2) THEN 
C RESERVED FOR ADS GRADIENTS 
GO TO 10 
EESE 
IF(INFO.EQ.-1) THEN 
C ADS OVERRIDE VALUES 
с (for no scaling, rel FD step, min |FD step|,zero) 
ТИК ( 2 )=0 


WK(21)=0.001D0 
WK (22 )=0.0001D0 
WK(26)=0.1D0 
WK(37)=1.0D-10 
GO TO 10 
ELSE 
C OPTIMIZATION COMPLETE, INFO-O 


IF(INFO.EQ.0) GO TO 20 
ENDIF 
ENDIF 
ENDIF 
20 CONTINUE 
C RECORD RESULTANT PRESSURE PROFILES 
DO 104 I=0,NTIME 
WRITE(89,*) TIME(I),PRESS(I),PCOMPAR(I) 
104 CONTINUE 
CLOSE(89) 
C PRINT DESIGN DATA TO SCREEN 
PRINT* 
PRINT*, ' FINAL DESIGN’ 
PRINT1001,(TAP(I),I-1,NTAP) 
1001 FORMAT(' LENGTHS = ЗЕ FT’) 
PRINT1002,(DIAM(I),I-0,NTAP) 
1002 ҒОВМАТ(! DIAMETERS = ',4F10.2,' IN’) 
PRINT1003 , WEIGHT 
1003 FORMAT(' CHARGE WEIGHT =',F6.1,' LB’) 
PRINT1005 , RANGE 
1005 FORMAT(' RANGE - SEE. qo one 
PRINT1004,0BJ 
1004 ҒОВМАТ(! SQRT OF AVG OF (Pdesign-Pdesired)^2 -', F6.1, 
c' PSI') 
PRINT* 
РЕТМТ1006 ,МСА 1,5 
1006 ҒОВМАТ(! CALLS TO P-T HISTORY GENERATOR =' ,I4) 
END 


П. SUBROUTINE DCOMPAR 


c SUBROUTINE DCOMPAR 

e This subroutine, for use with DTAPOPT, computes the square root of 
c of the average sum of the squares of pressure differences. 
с-------- 244 SUUM 2 „л = 
С PRECISION: DOUBLE 

c AUTHOR: ОА АТ rele le te eel 

c LAST UPDATE: 6/17/92 

c INTERFACE: 4 BUSES: OPTI, OPTDPA, FORDPTGEN, FORDCOMPAR. 


о о а – 
с VARIABLES 

C NAME ТҮРЕ COMMON OR LOCAL 170? DESCRIPTION 

C I I local na local index 

C NTIME I OPTI I calculated t 

C PCOMPAR(I,I=O.NTIME) DP FORDPTGEN I mud desired P 
C PRESS(I,I-2O,NTIME) DP OPTDPA I calculated P 

C SUMSQ pP FORDCOMPAR О value to be opt'd 
C Xe dee dede eese eee dee ee eee e eee dede eee dede eee eee de dede de ede e dee ee eee 


SUBROUTINE DCOMPAR 
IMPLICIT NONE 





C SPECIFICATIONS FOR GLOBAL VARIABLES 

COMMON/OPTI /NTAP , NTIME 
INTEGER NTAP , NTIME 
COMMON/OPTDPA/ ТАР, DIAM, TIME. PRESS 
DOUBLE PRECISION TAP(10),DIAM(0:10),TIME(0:1000) , PRESS(0:1000) 

| COMMON/FORDPTGEN/NDES , TDES , PDES, PCOMPAR 

| INTEGER NDES 
DOUBLE PRECISION TDES(0:1000) , PDES(0:1000) , PCOMPAR(0:1000) 


COMMON /FORDCOMPAR /SUMSQ 
DOUBLE PRECISION  SUMSQ 
С SPECIFICATIONS FOR LOCAL VARIABLE 
INTEGER I 
С BEGIN EXECUTION 
SUMSQ=0 . ODO 
DO 99 I=0,NTIME 
SUMSQ-SUMSQ- ( PCOMPAR( 1) -PRESS(I))*(PCOMPAR(I)-PRESS(I)) 
99 CONTINUE 
SUMSQ=SQRT( SUMSO/(NTIME+1 ) ) 
RETURN 
END 


ПІ. SUBROUTINE DPTGEN 


c SUBROUTINE DPTGEN 

C This subroutine, for use with DTAPOPT, interpolates the input 

с desired pressure-time history to find the desired pressure at each 

c time computed by the P-T generator. 

@ - = - - = - - - - = - - - - - - - - - - - - - - - - - - - - - - - - - 
ШРЕЕСІ5ТОМ: DOUBLE 

c AUTHOR: William Earl Miller II 

Guest UPDATE: 6/17/92 

с INTERFACE: о РОБЕ: OPTI, OFTDPA TRORDEPTCEN 

С = - - - - - - - - - - - = - - - = = = = - - - - - - = - - - - - - - - 
c VARIABLES | 

C NAME TYPE COMMON OR LOCAL 1/0? DESCRIPTION 

C I I local na local index 

С Ј | Тоса1 па local index 

e NDES I FORDPTGEN I No. input pts - 1l 
С NTIME I OPTI I No. calc pts - 1 
с PCOMPAR(I), I-O,NTIME DP FORDPTGEN О Int Desired Press 
С PDES(I), I2O,NDES DP FORDPTGEN I Input pressure 

С MDES(I), I=0,NDES DE FORDPTGEN I Input time 

C TIME(I), I=0,NTIME DP OPTDPA I Computed time 

C KKK KAKA KAKA AK AHH BAK KHAN AHA HANK AANA LE HAAN EAA KAHNE KK KKK KKH 


SUBROUTINE DPTGEN 
PHPLICIT NONE 
C SPECIFICATIONS FOR GLOBAL VARIABLES 


COMMON/OPTI/NTAP , NTIME 


INTEGER NIAE NTIME 


COMMON/OPTDPA/ TAP, DIAM, TIME. PRESS 
DOUBLE PRECISION TAP(10),DIAM(0:10), TIME(0:1000) , PRESS(0:1000) 
COMMON/FORDPTGEN/NDES , TDES , PDES, PCOMPAR 
INTEGER NDES 
DOUBLE PRECISION TDES(0:1000) , PDES(0:1000) , PCOMPAR(0:1000) 
SPECIFICATIONS FOR LOCAL VARIABLES 
INTEC HEROS J 
BEGIN EXECUTION 
1-0 


IF((ABS(TIME(J)).GE.1.0D-13).AND. (ABS(TDESCIOO СЕЛІ ОРБІР НИНЕ 
PRIN@, “INPUT PRESSURE PROFILE START LIME NOT Zuko 
STOP 
ELSE 
TIME(J)-TDES(J) 
ENDIF 
DO 100 I=0,NTIME 


1 ТЕСТІМЕСІ) СЕ IDES Ge then 


IF(TIME(I).LE.TDES(J+1)) THEN 
PCOMPAR(I)=PDES(J)+((TIME(I)-TDES(J))/(TDES(J+1)-TDES(J))) 
C *(PDES(J+1)-PDES(J)) 
GO TO 100 
ENDIF 
ENDIF 
J=J+1 
IF(J.EQ.NDES+1) THEN 
PRINTS" OUTISOESSENBESDSPEORSPOINBSE 
STOP 
ELSE 
GO TO 1 
ENDIF 


100 CONTINUE 


RETURN 
END 


IV. SUBROUTINE DTAPWT 


ОО ООО О Со О ООО О ОСОО О о-о оого 


SUBROUTINE DTAPWT 
This subroutine, for use with DTAPOPT, calculates the weight 
of a tapered charge. 


PRECISION: DOUBLE 


AUTHOR: William Earl Miller II 

LAST UPDATE: 6/17/92 

INTERFACE: 3 BUSES ОРТТ, ОРТБРА st co 

VARIABLES FROM COMMON STATEMENTS 
NAME TYPE COMMON OR LOCAL 1/0? DESCRIPTION 
DIAM(I),I-O,NTAP DP OPTDPA I Chg Diams (in) 

I I local na local index 
NTAP I ОРТТ I No. Chg Segs 

РІ DP local na Pi 

RHO DP local na H20 dens (lbm/ft^3) 
SGRAV DP local na Sp Grav of Chg 
TAP(I),I=1,NTAP DP OPTDPA I Length of Seg (ft) 
VOL DP local na Charge Vol 
WEIGHT DP PASS na Charge Wt (160) 


Vno eee eee eee eee ye RERRERRERRERRERERERERRERRERRERRRRERERRERRER 


SUBROUTINE DTAPWT 


2 





BMPLICIT NONE 


COMMON /OPTI/NTAP, 
INTEGER ҚТАР, 
COMMON /OPTDPA/ 
DOUBLE PRECISION 
COMMON/PASS/ 
DOUBLE PRECISION 


DOUBLE PRECISION 
INTEGER I 


SPECIFICATIONS FOR GLOBAL VARIABLES 
NTIME 
КТІМЕ 
ТАР, DIAM, TIME ERES S 
ПРИШАО МИ АЛА COMON TIME COT TOO0O0) ,PRESS(0:1000) 
WEIGHT 
WEIGHT 

SPECIFICATIONS FOR LOCAL VARIABLES 
EIT RHO, SGRAV, VOL 


BEGIN EXECUTION 


PI=4 .ODO*ATAN(1.0D0) 


RHO= 62.32D0 
SGRAV= 1.712D0 
VOL=0 . ODO 

DO 101 I-1,NTAP 


VOL=VOL+PI*TAP(1I)*(DIAM(I-1)*DIAM(I-1)+DIAM(I)*DIAM(I) 
c *DIAM(I-1)*DIAM(I))/1728.0D0 


101 CONTINUE 


WEIGHT=VOL*SGRAV*RHO 


RETURN 
END 
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